#ifndef _CU_FLOAT4_H_
#define _CU_FLOAT4_H_

//**********************************************************************************
// float4 operators and functions
//**********************************************************************************
SDH_CALL float4   operator +  (const float4 & u, float v) { float4 r = { u.x + v, u.y + v, u.z + v, u.w + v }; return r; }
SDH_CALL float4   operator -  (const float4 & u, float v) { float4 r = { u.x - v, u.y - v, u.z - v, u.w - v }; return r; }
SDH_CALL float4   operator *  (const float4 & u, float v) { float4 r = { u.x * v, u.y * v, u.z * v, u.w * v }; return r; }
SDH_CALL float4   operator /  (const float4 & u, float v) { float4 r = { u.x / v, u.y / v, u.z / v, u.w / v }; return r; }

SDH_CALL float4   operator +  (float v, const float4 & u) { float4 r = { v + u.x, v + u.y, v + u.z, v + u.w }; return r; }
SDH_CALL float4   operator -  (float v, const float4 & u) { float4 r = { v - u.x, v - u.y, v - u.z, v - u.w }; return r; }
SDH_CALL float4   operator *  (float v, const float4 & u) { float4 r = { v * u.x, v * u.y, v * u.z, v * u.w }; return r; }
SDH_CALL float4   operator /  (float v, const float4 & u) { float4 r = { v / u.x, v / u.y, v / u.z, v / u.w }; return r; }

SDH_CALL float4   operator +  (const float4 & u, const float4 & v) { float4 r = { u.x + v.x, u.y + v.y, u.z + v.z, u.w + v.w }; return r; }
SDH_CALL float4   operator -  (const float4 & u, const float4 & v) { float4 r = { u.x - v.x, u.y - v.y, u.z - v.z, u.w - v.w }; return r; }
SDH_CALL float4   operator *  (const float4 & u, const float4 & v) { float4 r = { u.x * v.x, u.y * v.y, u.z * v.z, u.w * v.w }; return r; }
SDH_CALL float4   operator /  (const float4 & u, const float4 & v) { float4 r = { u.x / v.x, u.y / v.y, u.z / v.z, u.w / v.w }; return r; }

SDH_CALL float4 & operator += (float4 & u, const float4 & v) { u.x += v.x; u.y += v.y; u.z += v.z; u.w += v.w; return u; }
SDH_CALL float4 & operator -= (float4 & u, const float4 & v) { u.x -= v.x; u.y -= v.y; u.z -= v.z; u.w -= v.w; return u; }
SDH_CALL float4 & operator *= (float4 & u, const float4 & v) { u.x *= v.x; u.y *= v.y; u.z *= v.z; u.w *= v.w; return u; }
SDH_CALL float4 & operator /= (float4 & u, const float4 & v) { u.x /= v.x; u.y /= v.y; u.z /= v.z; u.w /= v.w; return u; }

SDH_CALL float4 & operator += (float4 & u, const float v) { u.x += v; u.y += v; u.z += v; u.w += v; return u; }
SDH_CALL float4 & operator -= (float4 & u, const float v) { u.x -= v; u.y -= v; u.z -= v; u.w -= v; return u; }
SDH_CALL float4 & operator *= (float4 & u, const float v) { u.x *= v; u.y *= v; u.z *= v; u.w *= v; return u; }
SDH_CALL float4 & operator /= (float4 & u, const float v) { u.x /= v; u.y /= v; u.z /= v; u.w /= v; return u; }

SDH_CALL float4   operator - (const float4 & v) { float4 r = { -v.x, -v.y, -v.z, -v.w }; return r; }

SDH_CALL float3 xyz       (const float4 & v) { float3 u = { v.x, v.y, v.z }; return u; }

SDH_CALL float  dot       (const float4 & u, const float4 & v) { return (u.x*v.x + u.y*v.y + u.z*v.z + u.w*v.w); }
SDH_CALL float  dot       (const float4 & u, const float3 & v) { return (u.x*v.x + u.y*v.y + u.z*v.z + u.w); }
SDH_CALL float  dot       (const float3 & u, const float4 & v) { return (u.x*v.x + u.y*v.y + u.z*v.z + v.w); }
SDH_CALL float  dot3      (const float4 & u, const float4 & v) { return (u.x*v.x + u.y*v.y + u.z*v.z); }
SDH_CALL float  dot3      (const float4 & u, const float3 & v) { return (u.x*v.x + u.y*v.y + u.z*v.z); }
SDH_CALL float  dot3      (const float3 & u, const float4 & v) { return (u.x*v.x + u.y*v.y + u.z*v.z); }

SDH_CALL float4 pow       (const float4 & u, float g) { float4 r = { pow(u.x, g), pow(u.y, g), pow(u.z, g), pow(u.w, g) }; return r; }

SDH_CALL float4 clamp     (const float4 & u, float a, float b) 
{ 
  float4 r = { clamp(u.x, a, b), clamp(u.y, a, b), clamp(u.z, a, b), clamp(u.w, a, b) }; return r; 
}
SDH_CALL float4 saturate  (const float4 & u) 
{ 
  return clamp(u, 0.0f, 1.0f); 
}

SDH_CALL float  luminance (const float4 & v) { float3 u = {0.35f, 0.51f, 0.14f}; return (u.x*v.x + u.y*v.y + u.z*v.z); }

SDH_CALL float4 catmullrom (const float4 & P0, const float4 & P1, const float4 & P2, const float4 & P3, float t)
{
  const float ts = t * t;
  const float tc = t * ts;

  return (P0 * (-tc + 2.0f * ts - t) + P1 * (3.0f * tc - 5.0f * ts + 2.0f) + P2 * (-3.0f * tc + 4.0f * ts + t) + P3 * ( tc - ts)) * 0.5f;
}

SDH_CALL float4 hermitespline(const float4 & P0, const float4 & P1, const float4 & P2, const float4 & P3, float t, float c = 0.5f)
{
	const float ts = t * t;
  const float tc = t * ts;

	float4 m0, m1;

#define HERMITE_SPLINE_CARDINAL 1

#ifdef HERMITE_SPLINE_CATMULLROM
	m0 = (P2 - P0) * 0.5f;
	m1 = (P3 - P1) * 0.5f;
#endif

#ifdef HERMITE_SPLINE_CARDINAL
	m0 = (1.0f-c) * (P2 - P0) * 0.5f;
	m1 = (1.0f-c) * (P3 - P1) * 0.5f;
#endif

	return (P1 * (2.0f*tc - 3.0f*ts + 1.0f) + m0 * (tc - 2.0f*ts + t) + P2 * (-2.0f*tc + 3.0f*ts) + m1 * (tc - ts));
}

SDH_CALL float4 monotonebicubic(const float4 & P0, const float4 & P1, const float4 & P2, const float4 & P3, float t)
{
	const float ts = t * t;
  const float tc = t * ts;

	float dk0, dk1, dk2;

	float p0 = luminance(P0);
	float p1 = luminance(P1);
	float p2 = luminance(P2);
	float p3 = luminance(P3);

	dk0 = p1 - p0;
	dk1 = p2 - p1;
	dk2 = p3 - p2;

	float m0 = (dk1 + dk0) * 0.5f;
	float m1 = (dk1 + dk2) * 0.5f;

	float a = 0.0f, b = 0.0f, tau = 0.0f;
	
	if (abs(p1 - p2) > 0.25f)
	{
		a = m0 / dk1;
		b = m1 / dk1;

		tau = 3.0f / sqrt(a*a + b*b);

		m0 = dk1 * a * tau;
		m1 = dk1 * b * tau;
	}
	else
	{
		m0 = 0.0f;
		m1 = 0.0f;
	}

	float4 M0 = (P2 - P0) * 0.5f * m0;
	float4 M1 = (P3 - P1) * 0.5f * m1;

	return (P1 * (2.0f*tc - 3.0f*ts + 1.0f) + M0 * (tc - 2.0f*ts + t) + P2 * (-2.0f*tc + 3.0f*ts) + M1 * (tc - ts));
}

SDH_CALL float4 lerp      (const float4 & u, const float4 & v, float t ) 
{ 
  return u + t * (v - u); 
}

SDH_CALL float4 min (const float4 & u, const float4 & v) { float4 r = {min(u.x, v.x), min(u.y, v.y), min(u.z, v.z), min(u.w, v.w)}; return r; }
SDH_CALL float4 max (const float4 & u, const float4 & v) { float4 r = {max(u.x, v.x), max(u.y, v.y), max(u.z, v.z), max(u.w, v.w)}; return r; }
SDH_CALL float3 min3(const float4 & u, const float4 & v) { float3 r = {min(u.x, v.x), min(u.y, v.y), min(u.z, v.z)}; return r; }
SDH_CALL float3 max3(const float4 & u, const float4 & v) { float3 r = {max(u.x, v.x), max(u.y, v.y), max(u.z, v.z)}; return r; }
SDH_CALL float3 min3(const float4 & u, const float3 & v) { float3 r = {min(u.x, v.x), min(u.y, v.y), min(u.z, v.z)}; return r; }
SDH_CALL float3 max3(const float4 & u, const float3 & v) { float3 r = {max(u.x, v.x), max(u.y, v.y), max(u.z, v.z)}; return r; }
SDH_CALL float3 min3(const float3 & u, const float4 & v) { float3 r = {min(u.x, v.x), min(u.y, v.y), min(u.z, v.z)}; return r; }
SDH_CALL float3 max3(const float3 & u, const float4 & v) { float3 r = {max(u.x, v.x), max(u.y, v.y), max(u.z, v.z)}; return r; }

SDH_CALL float  sad (const float4 & u, const float4 & v) { return (abs(u.x - v.x) + abs(u.y - v.y) + abs(u.z - v.z) + abs(u.w - v.w)); } // Sum of Absolute Differences for xyzw components
SDH_CALL float  ssd (const float4 & u, const float4 & v) { return (sqr(u.x - v.x) + sqr(u.y - v.y) + sqr(u.z - v.z) + sqr(u.w - v.w)); } // Sum of Squared Differences for xyzw components
SDH_CALL float  sad3(const float4 & u, const float4 & v) { return (abs(u.x - v.x) + abs(u.y - v.y) + abs(u.z - v.z)); } // Sum of Absolute Differences for xyz components only
SDH_CALL float  ssd3(const float4 & u, const float4 & v) { return (sqr(u.x - v.x) + sqr(u.y - v.y) + sqr(u.z - v.z)); } // Sum of Squared  Differences for xyz components only

SDH_CALL float4 set(float4 & v, float x = 0.0f, float y = 0.0f, float z = 0.0f, float w = 0.0f) 
{ 
  v.x = x; v.y = y; v.z = z; v.w = w; return v; 
}
SDH_CALL float4 set(float4 & v, float3 xyz, float w = 0.0f) 
{ 
  v.x = xyz.x; v.y = xyz.y; v.z = xyz.z; v.w = w; return v; 
}
SDH_CALL float4 set(float4 & v, float2 xy, float z = 0.0f, float w = 0.0f) 
{
  v.x = xy.x; v.y = xy.y; v.z = z; v.w = w; return v; 
}

SDH_CALL float4 f4()                                   { float4 r = {  0,   0,   0,   0}; return r; }
SDH_CALL float4 f4(float x)                            { float4 r = {  x,   x,   x,   x}; return r; }
SDH_CALL float4 f4(float x, float y, float z, float w) { float4 r = {  x,   y,   z,   w}; return r; }
SDH_CALL float4 f4(const float3 & u, float w=0.0f)     { float4 r = {u.x, u.y, u.z,   w}; return r; }
SDH_CALL float4 f4(const float4 & u)                   { float4 r = {u.x, u.y, u.z, u.w}; return r; }
SDH_CALL float4 f4(const uchar4 & u)                   { float4 r = {u.x, u.y, u.z, u.w}; return r; }
SDH_CALL float4 f4_norm(const uchar4 & v)              { float4 r = { v.x / 255.0f, v.y / 255.0f, v.z / 255.0f, v.w / 255.0f }; return r;}

SDH_CALL float  length    (const float4 & u) { return sqrt(sqr(u.x) + sqr(u.y) + sqr(u.z) + sqr(u.w)); }
SDH_CALL float  length3   (const float4 & u) { return sqrt(sqr(u.x) + sqr(u.y) + sqr(u.z)); }
SDH_CALL float4 normalize (const float4 & v) { return v / length(v); }
SDH_CALL float3 normalize3(const float4 & v) 
{ 
	float3 u = {v.x, v.y, v.z};  
  float l = length3(v); 
  u.x /= l; u.y /= l; u.z /= l;
  return u;
}

SDH_CALL float4 normalize4(const float4 & v) 
{ 
	float3 u = {v.x, v.y, v.z};  
	float l = length3(v); 
	u.x /= l; u.y /= l; u.z /= l;
	return f4(u);
}

SDH_CALL float4 sqrt      (const float4 & u) { float4 r = { sqrt(u.x), sqrt(u.y), sqrt(u.z), sqrt(u.w) }; return r; }
SDH_CALL float4 sqr       (const float4 & u) { float4 r = { sqr (u.x), sqr (u.y), sqr (u.z), sqr (u.w) }; return r; }

#endif